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1 Introduction 

Frequency Map Analysis is a numerical method based on refined Fourier tech- 
niques which provides a clear representation of the global dynamics of many 
multi-dimensional systems, and which is particularly adapted for systems of 3- 
degrecs of freedom and more. This method relics heavily on the possibility of 
making accurate quasiperiodic approximations of of quasiperiodic signal given 
in a numerical way. In the present paper, we will describe the basis of the 
frequency analysis method, focussing on the quasi periodic approximation tech- 
niques. Application of these methods for the study of the global dynamics and 
chaotic diffusion of Hamiltonian systems and symplectic maps in different do- 
mains can be found in (Laskar, 1988, 1990, Laskar and Robutel, 1993, Robutel 
and Laskar, 2001, Robutel, 2003, this volume, and Nesvorny and Ferraz-Mello, 
1997) for solar system dynamics, and in (Papaphilippou and Laskar, 1996, 1998, 
Laskar, 2000, Wachlin and Ferraz-Mello, 1998, Valluri and Merritt, 1998, Mer- 
ritt and Valluri, 1999) for galactic dynamics. The method has been particularly 
successful for its application in particle accelerators (Dumas and Laskar, 1993, 
Laskar and Robin, 1996, Robin et ai, 2000, Comunian et ai, 2001, Papaphilip- 
pou and Zimmermann, 2002, Steier et ai, 2002), and was also used for the 
understanding of atomic physics (Milczewski et ai, 1997), or more general dy- 
namical system issues (Laskar et ai, 1992, Laskar, 1993, 1999, Chandre et ai, 
2001). 

2 Frequency Maps 

According to the KAM theorem (see Arnold et ai, 1988), in the phase space 
of a sufficiently close to integrable conservative system, many invariant tori 
will persist. Trajectories starting on one of these tori remain on it thereafter, 
executing quasiperiodic motion with a fixed frequency vector depending only on 
the torus. The family of tori is parameterized over a Cantor set of frequency 
vectors, while in the gaps of the Cantor set chaotic behavior can occur. The 
frequency analysis algorithm will numerically compute over a finite time span 
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a frequency vector for any initial condition. On the KAM tori, this frequency 
vector will be a very accurate approximation of the actual frequencies, while 
in the chaotic regions, the algorithm will still provide some determination of 
the frequency vector, but in this region, complementary of the KAM tori, the 
frequency vector will not be uniquely defined. 

Let us consider a n-DOF Hamiltonian system close to integrablc in the form 
H(I, 9) = HoW+eH^I, 9), where H is real analytic for (I, 9) e B n x T™, B n is 
a domain of M. n and T™ is the n-dimensional torus. For e = 0, the Hamiltonian 
reduces to H (I) and is integrable. The equations of motion are then for all 
j = l,...,n 

dHJI) 

ij = o . e i = -£rr L = 1 'jW ■> ( x ) 

which gives Zj(t) = Zjoe lVst in the complex variables zj = ijexpi#j, where 
zjo = Zj(0). The motion in phase space takes place on tori, products of true 
circles with radii Ij — |2j(0)|, which are described at constant velocity Vj{I). If 
the system is nondegenerate, that is if 

the frequency map F : B n — ► R"; (7) — ► (v) is a diffcomorphism on its 
image fi, and the tori are as well described by the action variables (I) G B n or 
in an equivalent manner by the frequency vector (y) G Q,. For a nondegenerate 
system, KAM theorem still asserts that for sufficiently small values of e , there 
exists a Cantor set f2 e of values of (v), satisfying a Diophantinc condition of the 
form 

\{k,y)\> n e /\k\ m (3) 

for which the perturbed system still possesses smooth invariant tori with linear 
flow (the KAM tori). Moreover, according to Poschel (1982), there exists a 
diffcomorphism 

* : TP x fi — ► T" x B n ; (p, v) — > (9, 1) (4) 

which is analytical with respect to <p, C°° in v, and on T" x Ct e transforms the 
Hamiltonian equations into the trivial system 

Vj = , (pj = uj . (5) 

For frequency vectors (v) in £ , the solution lies on a torus and is given in 
complex form by its Fourier series 

Zj (t) = zjoe^ + a m {vy <m - v>t (6) 

m 

where the coefficients a m {v) depend smoothly on the frequencies (y). If we fix 
9 G T" to some value 9 — 6q, we obtain a frequency map on B n defined as 

Fg :B n ^rt; I — ► p2( i &~ 1 (6o,I)) (7) 

where P2 is the projection on O (p2((/),v) = v). For sufficiently small e, the 
torsion condition (2) ensures that the frequency map Fg is a smooth diffcomor- 
phism. 
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2.1 Isoenergetic nondegeneracy 

If the isoenergetic nondegeneracy condition of Arnold (1989) is verified 

/ d 2 H (I) dH (I) \ 



det I 81 [ > 

V v{iy o 



det 



dF~ 

dH (I) ' 

V a/ 



dl 




then, if v n = — j— — ^ 0, and for small enough e, the application 



(8) 



{h, ■ ■ ■ ,I n -i,I n ) — > f — ; • ■ ■ 7 ,H 



(9) 



is a diffeomorphism, and so will be its restriction F' g on H 1 {K) 



(il, . . . ,l n -l) > , ■ • • , 



(10) 



When the isoenergetic nondegeneracy condition is verified, we can thus restrict 
ourself to an energy level H = h, and the frequency map Fg : M™ — > K" can be 
reduced to a map F' 0Q : R" _1 K n_1 . 



3 Quasiperiodic approximations 

The frequency analysis method and algorithms rely heavily on the observation 
that when a quasiperiodic function f(t) in the complex domain C is given nu- 
merically it is possible to recover a quasiperiodic approximation of f(t) in a 
very precise way over a finite time span [— T, T], several orders of magnitude 
more precisely than by simple Fourier analysis. Indeed, let 

/(*) = + J2 ake i{k ' v)i , a k e C (11) 
feez»-(i,o,...,o) 

be a KAM quasiperiodic solution of an Hamiltonian system in B n x T", where 
the frequency vector (y) satisfies a Diophantinc condition (3). The frequency 
analysis algorithm NAFF will provide an approximation f'(t) = J2k=i a'k^^ 
of f(t) from its numerical knowledge over a finite time span [— T, T] . The 
frequencies w' k and complex amplitudes a' k are computed through an iterative 
scheme. In order to determine the first frequency ui[, one searches for the 
maximum amplitude of (p(a) = {f{t) 1 e l<yt ) where the scalar product {f(t),g(t)) 
is defined by 

(f(t),g(t)) = ± J j{t)g{t)x{t)dt , (12) 

and where is a weight function (see next section). Once the first peri- 
odic term e 4 " 1 * is found, its complex amplitude is obtained by orthogonal 
projection, and the process is restarted on the remaining part of the function 
f 1 (t) = f(t)-a[e i ^ t . 
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In the next sections, we provide the rigorous foundations of the frequency anal- 
ysis algorithm by showing that this algorithm converges towards the frequency 
map described formally in the previous sections. It is interesting to see that 
the convergence of the frequency map algorithm will require some conditions on 
the frequency vector that will always be satisfied in the case of a KAM regular 
solution. 



4 Convergence of Frequency Map Analysis. 

Deflniton. 1 A weight function is a positive, even C°° junction \ on [ — 1 , 1] 
such that \ f_ 1 x(t)dt = 1. We will call the transform of x, the C°° function 
on R, Lp x , defined as 



<p x {x) = (e M , 1}? where (f(t), g{t))* = ^ J ^ f(t)g(t)x(t/T) dt 

(13) 



Lemma. 1 If x is a weight function, we have (p x (0) = 1; (p' x (0) = 0; f x {0) < 0. 

For all n G N, we have \(p x n \x)\ < 1, and there exists M n > such that 
\xip {n) {x)\ < M n onR. 

Proof. As x(t) is even, ty? x (x) is a real and even function, which implies tp x (0) = 
0. We have <p x (0) < as it is the integral of a strictly negative function. We have 

also \<pW(x)\ < \ j\ \t\ n X (t)dt < \ j\x(t)dt = 1. On the other hand, using 

integration by parts, we obtain |x<^™)(x)| < (x(l)+x( — 1))/2+5 J_ 1 x'{t)dt+n. 

Deflniton. 2 Let v = (1/1,1/2, ... ,v n ), and 

/(t)=e^*+ dke i{k ' v)t ; a k eC (14) 

fe€Z»-(l,0,...,0) 



a quasi periodic function onR. In the frequency map analysis, the approxima- 
tion i/f of Pi is obtained as the value of a for which f(o-) = | (f(t), e %at )\ | is 
maximum, in a neighborhood of V\ . 

Theorem. 1 Let \ be a weight function, and tp = ip x its transform with the 
asymptotic expressions when x — > 00 

(15) 

where n > 1 and where go (x) , gi (x) and g2(x) are bounded on K. Let f(t) be a 
quasi periodic function on the form (14), and for all k, fl k = (k,v) — v\; and 

assume that J2k IofI * s conver gent for p = 0, 1, and n. Then for T — > +00, 

vf — ► v\ and 

Ul ~ vI = ? W 9liQkT) + ° (t^) (16) 
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Proof. With x = (yi — cr)T, will be obtained for the maximum value of the 
modulus of 

ip(x)=ip(x)+ a k ip(x + n k T) (17) 

feez n -(i,o,...,o) 

or in an equivalent way, the maximum of its square tp(x)tp(x). This maximum 
thus fulfills the condition ip' (x)-ip (x) + ip{x) , 4> 1 (x) — 0, that is 



F{x, T) = <p{x)<p'{x) + ^(akai) <p(x + n k T)ip'(x + fi,T) 

k,l 

+ ^2 ft(a fc ) (<p(x)<p'(x + n k T) + tp'(x)tp(x + fi fc T)) = . 

k 

(18) 

Lemma. 2 Let x be a weight function, and ip = ip x its transform. Let us also 



assume that the series ^2 k \a k \ and J2 k 
Si. Then for all A > 0, 



are convergent with sum So and 



lim a k tp(x + fifcT) = 



(19) 



uniformly with respect to x 6 [—A, A]. 

Proof. Let J-(x, T) = J2 k a kf{x + fifcT) and e > 0. With a Taylor expansion of 
<p at order n } we have 



F(x, T) = y2 — Y] a k <pW(Sl k T) + — — - V a fe ^ (n+1) (&) where & 6 

(n + 1)! V 



p=0 P ' fc 

We have from lemma 1 

x 



„n+l 



n+1 



(n + 1) 



^afc^+^fc) 



A n+1 

-JnTTy So 



(20) 



(21) 



and this expression can be made arbitrarily small for sufficiently large values of 
n. On the other hand, with the notations of lemma 1 



i 

T 



^^(fifcT)^) (fifcT) 



< ^4 ■ (22) 



For any e > 0, there exists thus no £ N such that for all n > no, and all 

|^,T)| = -^— M p + e (23) 

and for T sufficiently large, the sum will be arbitrarily small, which ends the 
demonstration. 



Corollary. 1 With the hypothesis of lemma 2, we have 

lim F(x,T) = (p(x)ip'(x) 



(24) 



uniformly for x G [— A, +A\. The function F{x,T) can thus be continuously 
extended in a function F*(x, T) : [—A, A]x]0, +oo] — ► K. 
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We have also for the derivative 
dF(x,T) 



dx 



= <p(x)ip"(x) +f' 2 (x) 

+ K K) W{x)ip'{x + n k T) + (p(x)tp"(x + n k T) + tp"(x)tp(x + n k T)} 

k 

+ K(o fc ai) Mx + n k T)<p"(x + Q t T) + <p'{x + fl k T)<p'(x + fl t T)} . 

k,l 

(25) 



With the hypothesis of lemma 2, we will have 

dF(x, T) 



lim 



d:i 



<p(x)<p"(x) +f' 2 (x) 



(26) 



uniformly for x G [—A,+A]. The extension F*(x,T) thus has a continu- 
ous first derivative on [— A, A] x]0, +oo]. From Lemma 1 and (26), we have 
dF*/dx (0, +oo ) = f"(0) < 0, we can apply the implicit function theorem 
(Schwartz, 1992, p.176) to F*(x,T) at the point (0,+oo). Thus, there exists a 
neighborhood U of +oo, and a unique continuous map x{T) = x such that 



^lim x(T) = 0; F(x(T), T) = for all T <= U . 



(27) 



Once the existence of x(T) verifying (27) is known, we can obtain a generalized 
equivalent for x(T) when T — > +oo. 

Lemma. 3 Let ip(x) be a function on K ; and n G N,n > 7 and assume that 
for all < p < n there exists M p > such that \x p ip(x)\ < M p for all x G K. 
TTien, /or a// A > 7 and all < p < n, there exists N p > suc/i t/iat 



|X*V(^ + ^)I<^ P V(x,X) G [-A, +A] 



(28) 



Froo/. We have + X) = (X + x) p (p(x + X) - ££0 C k X k x p ^ k ip(x + X) 

and the proof is obtained by recurrence on p. 



Lemma. 4 Let ip(x) be a C°° function on R, and assume that 



<p(x) 



9o(x) 



(29) 



where go(x) and gi(x) are bounded on R, and that the series J2 k 
vergent. Then for T — ► +00, 

J] a t ^(T) + n fc T) = J ^ft.(n»T) + [ i j . 

k k k \ / 



Proof. Let 



is con- 



(30) 



^(T) = T™ ( £ ak¥>{x(T) + n k T) ]T _^_ fl0 (n fcT )) ) . (31) 

\ fe k k / 
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Thus 

k k 

with T k {T) = (fl k T) n ip(x(T) + n k T) - g (n k T). From lemma 3, |^ fc (T)| is 
bounded, and the series (32) is uniformly convergent with respect to T. We 
thus have 

Urn HT)=Y / ^ lim T k {T) (33) 

k K 

But <p(x(T) + r> fc T) - ^(Q fe T) = x{T)y'{£ k {T)) with |&(T) - VL k T\ < \x(T)\. 
Thus, as x n (p'(x) is bounded on R, we have for all k 

lim T k {T) = 

1™ ^^T^OT^CO) + (fi fc T)>(n fc T) - go (n k T) = , 

which ends the proof of the lemma. We can now complete the proof of theorem 
1. Assuming the hypothesis of theorem 1, we can apply this lemma to ip and 
<p'. We have thus 

a k <p(x(T) + n k T) = J2 ^9o(n k T) + o(^) 

k k k \ / 

Y j a k p\x{T) + n k T) = Y J ^gi^kT) + ol^\ . 

k k k V / 

We have also ip(x) = l + o(x), and ip'(x) = ip"(0)x + o(x), and limT^ +00 x(T) = 
0. We can thus make the expansion of the expressions involved in Eq.18 and 
obtain, by an expansion at first order in x and order n in 1/T 

p"(0)x(T) + o(x(T)) + ^9i(n k T) + o = , (36) 

from which the proof of the theorem ends easily. Indeed, 

T n x{T) = -L^ (Y ^Wt) + 0(1)) , (37) 



(35) 



¥>"(0) + 



x(T) 



As lim T ^ +oc z(T) = 0, and <p"(0) ^ 0, then T n x(T) is bounded and thus 
o(x(T)) = o(l/T n ) from which 



*(T) 



?^)n»7™ fll(nfcT) + (^) ' (38) 



This end the proof of theorem 1. We can then prove the following result 
Theorem. 2 With the hypothesis of theorem 1, for T — ► +oo 

v^-vl = Y RKV + 0(1/^+1) (39) 
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This improvement requires the following lemma : 

Lemma. 5 Let \ be a weight function, ip — ip x its transform, and assume that 
< ^ 9°( X K ( X \ >r ^ gi( - x K I 1 ^ " ( , 92(x) ( 1 \ 



where go(x), gi(x), 32(2;) are bounded on R, and that the series J2k 
vergent for p = 0,n. Then for T — > +00, 

J2 a k <p(x(T) + n k T) = awtflkT) + 0{l/T 2n ) 
k k 

]T a W '(x(T) + n k T) = a k <p'(n k T) + 0(l/T 2n ) 



(40) 
is con- 



(41) 



Proof. We will make the proof of the first relation (for f(x)); the argument for 
f'(x) is the same. Let 

x(T) 2 

ip(x(T) + n k T) = v (n k T) + x(T)v'(n k T) + ^-L^(( k ) (42) 

where Cfc £ [^ k T, Sl k T + x(T)], and 

HT) = T 2n a k Mx(T) + n k T) - V (Q k T)]\ (43) 



that is, 



T(T)=T 2n J2a k 



x(T) 2 

x {T)y'{Sl k T) + ^J-v"{Qk) 



E^CH ( 44 ) 



From theorem 1, \T n x{T)\ is strictly bounded by a constant N. Moreover, from 
the hypothesis of the lemma, x n ip'(x) is also bounded by a constant N', and 
thus \n^T n ip'(n k T)\ < N', and 



\a k T 2n x(T)<p'(n k T)\ < NN' 



a k 



(45) 



Finally, \<p" (x)\ < 1 (lemma 1), and as \T 2n x(T) 2 \ < N 2 



J- a k — - — if (Cfc) 



N 2 
< — \ak\ 



(46) 



In (44) \T k {T)\ is thus bounded by u k = NN' 
general term of a convergent series of sum S, and 



a k 



N 

\a k \ which is the 



\HT)\<Y,\Fk(T)\<J2u k = S 



(47) 
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which ends the proof of the lemma. The proof of the theorem can now be 
completed. 

Proof of theorem 2 : From theorem 1, we already know that when T — > +00, 
the solution x(T) of EQ. (18) tends to zero and from (38), 

x(T) = 0(1/T n ) (48) 

We have thus also the expansions for T — > +00 

<p{x(T)) = l + 0(l/T 2n ) 



<p'{x{T)) = <p"(0)x(T) + 0(1/T 3n ) = 0(l/T n ); 
while Eq.(35) gives 

^^(1(11+^ = 0(1^) 

k 

J2aw'(x(T) + n k T) = 0(l/T n ) . 



(49) 



(50) 



Using these expressions, we can now expand all expressions in Eq.18 up to 
0(1/T 2n ) which gives 

v"(0)x(T) + ]T K{a k W{x{T) + Q k T) = 0(1/T 2n ) (51) 
k 

and with lemma 5 

x(T) = X{a k )<p'{il k T) + 0(1/T 2 ") (52) 

V (°) k 



4.1 Amplitudes 

Once the estimate of the precision on the frequencies is obtained, the precision 
on hte amplitudes is easily calculated. Indeed, with the same notations, the 
amplitude of the first term of f{t) in eq.14 (which exact value is 1), will be 

a\ =< f(t), c<* >= <fi(x(T)) + a k cp(x(T) + n k T) . (53) 

k 

Using the above estimates, one then finds that the error on the amplitude is 

a i - 1 = E a Mn k T) + 0(1/T 2n ) = 0(1/T") . (54) 
k 



5 Cosine windows 

We can now apply this theorem to specific examples of weight functions. The 
most simple weight function will be defined on [—1,1] as xo(t) = E with the 
associate transform 

sin a; 

( Pxo( x ) = (55) 
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-5 5 10 

Figure 1: Graph of <fo{%) and <Pi(x). 



With this window (or absence of window) , the decreases of the sidelobes of the 
transform is slow, (as l/x). This is due to the discontinuity at —1 and +1 of the 

window function Xo( x ) = lr i,i] 5 considered as a function on R. In this case, the 

transform <p X o( x ) * s the usual Fourier transform of xo{x),and it is known that 
the decrease at infinity of tp Xo ix) depends on the class of regularity of Xo(t). 
More precisely, 

Proposition. 1 If f{t) is of class C p on R, absolutely integrable, with all 
derivatives absolutely integrable, then its Fourier transform f(x) satisfies 



/ + OC 
-OC 



f{t) e ixt dt = oil/x p ) 



when 



+oo 



(56) 



We will thus improve the decreasing at infinity of the transform function ip x (x) 
by increasing the class of regularity of the extension to R of the window function 
%(£). This can be done by using the window function Xi(*) = l + cos7ri. Indeed, 
we have then x'i(*) = — 7rsin7ri, xi' W = 71-2 cos7ri and thus Xi(l) = Xi(— 1) = 
Xi(l) - Xi(-l) = 0, while x'i'(l) - xi'(-l) = The extension xi(t) is thus 
C 1 and the decreasing at infinity of the sidelobes of the transform <p Xl (x) will 
be much stronger (Fig.l). More generally, we can search for a trigonometric 
polynomial Xv sucn that the extension x P (i) is of class C p over R. Indeed, 



Proposition. 2 For all p e N, let x P (t) = n^r{ 1 + cos irt) p . Then x P {t) 
is the unique trigonometric polynomial P(cos irt, sin7rf) of degree < p that is a 
weight function of class C 2p ~ 1 . Its associated transform is 



ip P ix) 



- f 



e lxt x P (t) dt 



i-l) p 7T 2p ipl) 2 s\ 



I 1 2 
X(X — 7T 



)-(^-pV) 



(57) 



Proof. As a weight function is even, P can be expressed uniquely as a polynomial 
in cos7rf. Setting u(t) = cos7rt, the derivatives of P(«(t)) are given by the Faa 
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di Bruno formula 

p(u) {k) = t E ^k^ p{r){u) - uikl) -- u{kr) (58) 

r=l fcH \rk r =k ' 1 ' T " 

where fci, . . . , k r > 1. As for fc > 0, 

u (2k ){t) =7r 2 fe( _ 1)fecos(7rf)) 

u (2fc+i)( t ) =7 r 2fe+1 (-l) fe+1 sin(7rf); 

For fc > 0, we have u( 2fc+1 )(l) = u^ k+1 ^(-l) = 0, and u^(l) = u^(-l) = 
(— l) k+1 . The only index fci for which the contribution in Eq.58 is not zero are 
thus even, and we obtain 

P( u )( 2fe+1 )(l) = = 

(60) 

To Assume that the weight function is of class C 2p 1 is equivalent to assume 
that P(u)W(l) = P(u)W(-l) = for fc = 0, . . . 2p - 1, that is, from (60), 
p( fe )(— 1) = 0, for fc = 0, 1, . . . ,p — 1. P is thus a polynomial of the form 
P{X) = (1 + X) P Q(X), where Q(X) is a polynomial, and the proposition 
follows. The constant is determined such that 1/2 ^ 1 x p (t)dt = 1. We have 
also easily the asymptotic expansions for x — > oo 

(-ljV^pQW"^ / 1 \ 

< } (*) = + °l^+r) ; ( 61 ) 



and the expansion at the origin 



{x ) = l + ^x>+o(x 2 ) with 4'(0) = --^£-f^ (62) 



Thus, if /(i) is a quasi periodic function of the form (14), for which J^ fe 
is convergent for m = 0, 1, and 2p + 1; for T — > +oo, 



1/1 ~ 1/1 = ^(0)T 2p+2 E ftp 1 C0S ^ fcT ) + ° 1^+2 I ( 63 ) 
Using theorem 2, we have as well for the cosine windows ip p (x), 

Vl ~ ^ = IVW) ^ R (°*K( fi * T ) + °(Vr 4p+2 ) (64) 

The computation of the derivative y p (a; ) is easy if we write v?p(:e) = Csina;/D(a;) 
where 

C = (-l) p 7r 2p (p!) 2 ; D(a;) = x(a; 2 - tt 2 ) • • • (x 2 -pV) (65) 



11 



Then 



with 



<p'(x) 



C 
D{x) 



cos a; — sin a;- 



D(x) 



2x 



1 

^ — 2 2 
X X — 7T 



+ ••• + 



2x 

— 2 S - ? 

X — p 7T 



(66) 
(67) 



5.1 Exponential window 

The cosine weight function of order p, Xp(t) provides a window function of 
class C p , and thus improves the convergence of the frequency map algorithm 
(theorem 1 and 2). Although it will not be always very useful in practice, we can 
also consider C°° windows functions, for which x^ n \~ 1) = X^iX) = f° r au n > 
as for example x*(t) = cexp(-l/(l - t 2 )), with c w 0.22199690808403971891. 
In this case, we have for all values of n,m, lim^+oo a;™^ m ^(a;) = 0. Thus, if 
J2k In™ I are convergent for all m, as for a KAM solution (see below), we will 

k 

have for all n, limT^+oo T n {v\ — v{ ) = 0. Let us note that for this exponential 
window, MO) « -0.035100738376487704994. 



5.2 KAM solutions 

If /(t) is a KAM solution, with rotation vector (i/), that fulfills a a diophantine 
condition (3), it can be expressed as a quasiperiodic series of the form (6), 



f(t)^Y, a ^ <kfi> with 9i = Vi t 



(68) 



that is analytical with respect to the angles (6). The series ^2 k are then 



convergent for all p. Indeed, we have then 



a k 



\a k 



<^(l + |fc|) pm - 



(69) 



But due to the analyticity of the series (68) with respect to 8, there exists a 
small s > such that 



I a/el < 



(70) 



As \k\ = | fci| + • • • + \k n \, where n is the number of degrees of freedom of the 
Hamiltonian system, 



(1 + 1*1) < (l + N)(l + N)...(l + |fcn| 



E W ^ E ^(1 + |fcil) pm ■••(! + |fc«D Pm e- 

fel,...,fe„ Ke 



|fc|s 



(71) 



(72) 



= i T i n E( i+ i^i) pmeHfcsis - 

K e i=l,...,nfe i >0 

These latest series are convergent which ends the proof. In fact, this case will be 
the most useful, as the solutions we are looking for will usually be the quasiperi- 
odic KAM solutions of an analytical Hamiltonian system. In this case, using 
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the weight \ p , we can obtain an accuracy for the frequency analysis of order 
1 /T 2p+2 on these solutions, which gives in particular 1 /T 4 when using the Han- 
ning window (xi). 

6 Quasiperiodic decomposition 

The previous results tell us that with the proposed algorithm, we will recover 
in a very efficient manner the frequency v\ of the leading periodic term of the 
quasiperiodic decomposition of 

k 

A first approximation of the corresponding amplitude A\ will be given by 

4 1} =< /,ei > where e x = e lVlt (74) 

we are then left with the remainder 

/i = /-</, ei >ei (75) 

and we can start the process again in order to search for the frequency v-i of the 
next term. A technical complexity appears here, as e\ = e Wl * and e 2 = e W2 ' 
are not orthogonal for A classical way to overcome this difficulty is 

given by Gramm-Schmidt orthogonalization process. Indeed, if (ei, . . . ,e„) is 
an independent family of unit vectors, one can construct an orthogonal family 
(e[, . . . , e' n ), such that for all 1 < k < n, Vec(e\, . . . , e k ) = Vec(e[, . . . , e' fe ), 
where Vec(ui, . . . , Uk) denotes the vector space generated by (tti, . . . , Ufe). In- 
deed, 

n-l 

e n = e«/ ||e„|| where e„ = e n - < e„, e' k > e'k , (76) 

fc=i 

and 

fn = fn-l- < fn-l,e' n > e' n . (77) 

As f n -i _L Vec(ei, . . . , e„_i) = Vec(e[, . . . , e^j), it is interesting to note, for 
practical computation that 

< f„-i,e' n >=< / n _i, e„ > . (78) 



7 Numerical simulations 

The asymptotic convergence of (vi — v"[) when T — > +00 provides a good 
indication of the possibilities of the method, but in practice, this asymptotic 
behavior will be also limited by the value of the involved constants, and by 
numerical accuracy. In the following, the previous asymptotic expressions are 
tested for a very simple example of a quasiperiodic function. Instead of using 
the output of a numerical integration, we prefer here to take directly a quasi 
periodic function with 2 independent frequencies. The tested function is 

1 

*i(*) = =7 -r (79) 

W l + l/2c lt + l/4c-^* V ; 
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with lu = 2.02. This function is chosen because of the slow convergence of its 
Fourier expansion. The results are displayed in figures 2-6 and Table 1 with a 
cosine window Xp of various order p = 0, 1, 2, 3, 4, 5, and an exponential window. 
A second test function F 2 (t) is also constructed with the first 50 periodic terms 
of the frequency decomposition of F\ (t) . 



7.1 Corrections 

Theorems 1 and 2 provide equivalents of the frequency error (v\ — v[) when 
T — ► +oo. If a quasiperiodic approximation of f{t) is known, for example by 
a first frequency analysis, then these estimates can be used as corrections for 
frequency. It should be noted that when the correction is searched, one does not 
know the exact values of the perturbing frequencies involved in the formulas 
of theorem 1 and 2, but an approximate value Q' k . Under the hypothesis of 
theorem 1, one has 

n k = n k + o(^^J (so) 

and the expression (16) is as well valid with £l' k instead of Cl k . In the same way, 

(81) 



v '(n' k T) = cp'(n k T) + o(^j (82) 

and the estimation of theorem 2 is also valid with Vt' k instead of Cl k . This allows 
to simplify greatly the use of these corrections as the first estimation is sufficient 
in general. 




7.2 Discussion 

We have first analyzed the results of the frequency analysis of Fi (t) for differ- 
ent windows Xp> an d the results of the precision of the determination of the 
frequency versus the length of the time interval T are given in Figure 2. For 
each case, the results are fitted with a line and the the results are gathered in 
Table 1 where ao (column 5) is the slope of the straight line least squares fitted 
to the curves of Fig. 2. The same study is then repeated with a correction step 
provided by theorem 1 and formula (63) using 10 and 50 terms . The results 
are given in Figures 3-4 and in column a w and 050 of Table 1. These values 
need to be compared with the theoretical value obtained with the full correc- 
tion of theorem 1, given in column a^. The correction obtained with theorem 
2 (formula (64)) is then given in column a' 50 , with a correction computed also 
with 50 terms. 

It should be noted that the correction step improves the determination of the 
frequency, but only to a certain limit. With F\(t), it is usually hopeless to expect 
the approximation in 0(1/T 2n+1 ) given in theorem 2 (column a'^ of Table 1). 
Indeed, because of the slow decreases of the amplitude of the quasiperiodic 
expansion of Fi(t), the contribution of the neglected terms will dominate beyond 
a certain level of precision. 
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We have then performed the same experience with F%(t). In this case, the cor- 
rections performed with 50 terms can be considered as optimal, as the correction 
step is now performed with a very precise estimate of the remaining terms. In- 
deed, the results that are plotted in Figs 6,7, and gathered in columns bo, 650, 650 
of Table 1 are now in much better agreement respectively with the theoretical 
coefficients a, a^, a'^. 

It is clear from these numerical simulations that even for a perfectly quasiperi- 
odic function, there are some limitation on practical computation of the fre- 
quencies, especially when the quasiperiodic expansion is slowly converging. It 
is nevertheless striking to see that even in the case of a slow convergence, the 
agreement of a and ao is excellent. When the correction is used, there is still a 
good agreement although not as good, as the effect of the neglected terms will 
dominate after a certain level of precision, unless, as with F^[t), the quasiperi- 
odic expansion of the considered function is fully recovered with good accuracy. 



p 


—a 




~ a 'oo 




-aio 


— a 50 


~ a 50 


-b 




-&50 





2 


3 


3 


1.98 


2.11 


2.97 


2.94 


1.91 


3.15 


3.11 


1 


4 


5 


7 


3.94 


4.25 


4.90 


4.43 


3.98 


4.81 


6.04 


2 


6 


7 


11 


5.82 


6.45 


6.35 


6.11 


5.89 


7.01 


12.68 


3 


8 


9 


15 


7.64 


7.99 


8.01 


7.87 


8.69 


9.49 


18.85 


4 


10 


11 


19 


8.89 


9.01 


9.03 


8.71 


11.42 


12.48 


23.49 


5 


12 


13 


23 


9.15 


9.33 


9.33 


8.87 


13.50 


14.47 


26.24 


* 








5.46 








5.43 







Table 1: Slope of the fitted errors on the frequency, p is the order of the window, 
a the exponent of the error from eq. (63), ao the fitted slope with no correction, 
aio and 050 the correction using 10 or 50 periodic terms and the theoretical 
slope when all terms of equation (63)are used for the correction (providing the 
series is converging) for the determination of the first frequency of F\(t) (eq.79). 
a'^ is the same with equation (64). &o>^50)^5o are the corresponding value for 
the function F 2 (t) obtained as a truncated expansion of F\ (t) with 50 periodic 
terms. 

Overall, it appears that when the quasiperiodic expansion is not fully recovered 
(up to a certain precision), the correction obtained through equations (63) and 
(64) does not improve very much the results that can be deduced from the use 
of a window of higher order (Table 1). Considering that the computation of 
many terms is time consuming, and that it is not so easy to obtain always the 
same number of relevant terms for various initial conditions corresponding to 
different regularity of the solution, I would rather recommend to use different 
orders for the window of the data. The more the solutions are regular, the 
higher one can choose the order of the window. Practically, this order is still 
limited to p = 3, 5 for quasi periodic function, and maybe lower for numerical 
solutions of dynamical systems. When the motion is more chaotic, p = 1, or 
even sometimes p = may be preferred. In practice, one would increase the 
order of the window until the precision seems to decreases and use the highest 
possible value for p. 
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7.3 Exponential window 

We have also tested the exponential window of section (5.1), and the results 
are reported for p = * in Table 1, but as it can be seen in Table 1, this was 
not very successful, compared to the cosine windows. Investigation of different 
exponential windows should probably be continued, as it is not clear why such 
windows could not attain the accuracy reached by the cosine windows \p- 
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Figure 2: Error on the measured frequency versus duration of the integration 
for the first frequency of Fi(t), log 10 (err) is plotted versus log 10 (T/27r). x p is 
the result for a window of order p . 
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Figure 3: Error on the measured frequency versus duration of the integration 
for the first frequency of Fi(t). log 10 (err) is plotted versus log 10 (T/27r). x p is 
the result for a window of order p with a correction using equation (63) with 10 
terms. 
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Figure 4: Error on the measured frequency versus duration of the integration 
for the first frequency of Fi(t). log 10 (err) is plotted versus log 10 (T/27r). x p is 
the result for a window of order p with a correction using equation (63) with 50 
terms. 
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Figure 5: Error on the measured frequency versus duration of the integration 
for the first frequency of Fi(t), log 10 (err) is plotted versus log 10 (T/27r). x p is 
the result for a window of order p with a correction using equation (64) with 50 
terms. 
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Figure 6: i*2(i)- Error on the measured frequency versus duration of the inte- 
gration for the first frequency of F 2 (t). log 10 (err) is plotted versus log 10 (T/27r). 
x p is the result for a window of order p with a correction using equation (63) 
with 50 terms. 
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Figure 7: F2(t). Error on the measured frequency versus duration of the inte- 
gration. log 10 (err) is plotted versus log 10 (T/27r) for the first frequency of F 2 (t). 
x p is the result for a window of order p with a correction using equation (64) 
with 50 terms. 
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8 Beyond Nyquist frequency 
8.1 Nyquist frequency 

A typical limitation in spectral analysis is given by the so called Nyquist fre- 
quency. Roughly speaking, it will not be possible to determine a frequency v 
larger than n/h, where h is the sampling time interval of the observed data. If 
| | > n/h, the observed frequency will become v = v>q + kin /h, where k G Z 
such that \v\ < n/h. This is actually what can be observed using the NAFF 
algorithm detailed above. 



z/o/tt 


v/n 


(v - V)/-K 


0.990 


0.990 


0.000000000000000000 


0.991 


0.991 


0.000000000000000000 


0.992 


0.992 


0.000000000000000000 


0.993 


0.993 


0.000000000000000000 


0.994 


0.994 


0.000000000000000000 


0.995 


0.995 


0.000000000000000000 


0.996 


0.996 


0.000000000000000000 


0.997 


0.997 


0.000000000000000000 


0.998 


0.998 


0.000000000000000000 


0.999 


0.999 


0.000000000000000000 


1.000 


-1.000 


2.000000000000000000 


1.001 


-0.999 


2.000000000000000000 


1.002 


-0.998 


2.000000000000000000 


1.003 


-0.997 


2.000000000000000000 


1.004 


-0.996 


2.000000000000000000 


1.005 


-0.995 


2.000000000000000000 


1.006 


-0.994 


2.000000000000000000 


1.007 


-0.993 


2.000000000000000000 


1.008 


-0.992 


2.000000000000000000 


1.009 


-0.991 


2.000000000000000000 



Table 2: The function exp(zi'oi) is evaluated for t = —1000, 1000, with a stepsize 
h = 1. For fo/ir < h, the recovered frequency (V) is recovered by NAFF up 
to the machine precision, but for v§/n > h, the recovered frequency is here 
z/q — 2n/h. 



8.2 Multiple timescales problem 

The Nyquist aliasing constraint means that to recover a given period, one needs 
to sample the data with at least two points per period. On the opposite, in order 
to determine precisely the long periods, one needs that the total interval length 
T is several time larger than these periods, in order to reach the asymptotic 
rates of theorems 1 and 2, or to be able to separate properly close frequencies. 
One thus realizes that if good low frequency determination imposes that T is 
large, while high frequencies impose that h is small, we will face a problem when 
two very different time scales are present in the system. This is actually the case 
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for planetary systems where short periods (of the order of the year) are present, 
as well as long secular periods, ranging, in the Solar system for example, from 
40 000 years to a few millions of years (Laskar, 1990). 

The first frequency analysis of the solar system solutions were made uniquely 
on the output of the secular equations, where only the long periods were present 
(Laskar, 1988), but if we want now to perform a frequency analysis of the direct 
output of a numerical integration of Newton's equations, without filtering or 
averaging, both time scales will be present. 

For the Solar system, for example, the determination of the long secular frequen- 
cies with a good accuracy will require that T is larger than 20 millions of years, 
while the sampling h must be smaller than half of the shortest period of the 
system, that is a few days if we consider the Moon motion, while h = 5000 years 
was enough for the secular equations. The amount of data to handle becomes 
then considerable as well as the related numerical computations. 



8.3 Multiple timescales solution 

The solution that can be use to overcome this problem is simply to sample the 
data with two different sampling intervals that are very close h and h! = h + s. 
For a real x, we will denote [x] the integer such that 

- 1 -<x-[x]<y (83) 

With this notation, the frequency analysis of f(t) = cxp(i2nv t) will give a 
frequency 2nv such that 

uoh = vh + k with k = [voh] . (84) 

The Nyquist condition is then expressed by the fact that \v$h\ < 1/2 implies 
fc = 0. We assume now that the second sampling time interval h' = h + e is 
such that 

\v s\ < - . (85) 



We have then 

with 

and thus 

As \v Q s\ < 1/2, 

and thus, from (88), 



v h' = v (h + e) = vh + ev a + k 
= v'{h + e) + k' 



(86) 



k' = k + [vh + evq] (87) 
v'h! — vh = ev — [vh + ev ] . (88) 

[vh! - vh] = -[vh + ev ] , (89) 



v = {v'h 1 -vh - [v'h 1 -vh]) , (90) 

h - h 
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which allows to revover in all cases the true value of the frequencies. Practically, 
although this formula can be used, there is an easy way to reduce the numerical 
errors of these computations. Indeed, it is easy to see that 

k=— h ((„> - U )h' -[v'h' -vh]) , (91) 
h — h 

which allows to recover vq through (84) . The main advantage of this intermedi- 
ary step is that as k should be an integer, the computed value can be rounded 
to its nearest integer and 

\k] 

vo = v + V • (92) 



8.4 Numerical examples 

8.4.1 Single periodic term 

In order to evaluate numerically the efficiency of the method described above, we 
first used a function f(t) — exp(ifo t) with a single periodic term of frequency 
vq. f(t) is evaluated from —1000 to +1000, with a stepsize h = 1, and the 
frequency analysis is performed with different values of vq, for values of vq/tt 
higher than 1000. The value vq is only recovered for vq/it < 1 (Table 2), but for 
larger values of v , the use of a second stepsize hi — 1.001 and formula (91,92) 
allows to recover the true frequency for values of vo/ir as high as 1000 with great 
accuracy (Table 3). 

8.4.2 Sun-Jupiter-Saturn system 

As a second, and more realistic example, we have consider a full numerical 
integration of the Sun- Jupiter-Saturn system over 50 millions of years. The 
numerical integration is performed with the symplectic integrator SBAB^ that 
is well adapted to perturbed Hamiltonian systems (Laskar and Robutel, 2001). 
In order to illustrate the previous section, we have performed the frequency 
analysis of the variable z§ — exp(iw5) (Fig. 4), where e$ is the eccentricity of 
Jupiter, and ©5 its longitude of perihelion, and the analysis of sin(is/2) enp(iflz) 
(Fig. 5), where i$ and are the inclination and longitude of the node with 
respect to the ecliptic and equinox J2000 reference frame. If one considers the 
invariance of the angular momentum, this problems has 5 degrees of freedom, 
that should correspond to 5 fundamental frequencies for a regular KAM solution. 
These frequencies will be n^,n%, the mean mean motion of Jupiter and Saturn, 
55) 96 related to the precessional motion of the perihelion of Jupiter and Saturn, 
and sq related to the motion of the node of the two orbits (see Laskar, 1990). In 
the integration, we used an output stepsize of 200 years, while the integration 
stepsize is 0.1 year. This allows to recover precisely the secular frequencies of 
the system, but as the short period perturbations are important, many terms 
appear in the quasiperiodic decomposition that are in fact aliased terms coming 
from the short period terms (see Table 4,5). 

When we use an additional output time hi = 200.0002 years, and formula 
(91,92), we can recover the true value of the aliased frequencies voi, even if 
some of the periods are smaller than 6 years (columns Pi). In Table 4, 5, k is 
the integer appearing in formula (84), denoting the number of turns that have 
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Table 3: The function exp(ifoi) is evaluated for t — —1000, 1000, with a stepsize 
h = 1, and with bl = 1.001. The reconstruction formula (91,92) allows to obtain 
the correct frequency up to v^jix = 1000. 
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Table 4: Frequency analysis of 05 = es expfWs) in the Sun- Jupiter-Saturn 
system. The integration of the complete Newton equations is performed over 
50 myr with two output stepsizes h = 200 yr and h! = 200.0002 yr. v oi are 
the reconstructed frequencies using formula (91,92) and the two step sizes h 
and h! . Pi is the period of the terms, while Ai their amplitude. The different 
frequencies are identified as integer combination s of the fundamental frequencies 
voi = fcii«5 + k 2l n 6 + k 3l g 5 + k^ge + k 5i s 6 
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Table 5: Same as Table 4 for sin(z5/2) exp(ifi5), where i$ and fis are the incli- 
nation and longitude of the node with respect to the ecliptic and equinox J2000 
reference frame. The constant term is due to the invariance of the angular mo- 
mentum. It would not be present in the reference frame of the invariant plane, 
orthogonal to the angular momentum of the system. 
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Table 6: Fundamental frequencies of the Sun- Jupiter-Saturn system, obtained 
by frequency analysis over 50 myr. 

been "lost" by the large stepsize. All the determined frequencies v oi can then 
be identified as integer combinations 

v^i = k u n 5 + k 2i n 6 + k 3i g 5 + huge + k 5i s 6 (93) 

of the fundamental frequencies 715, n^, 35, s 6 given in Table 5. The full solu- 
tion can then be compared to the quasiperiodic solution obtained by iteration 
of (Bretagnon and Simon, 1990). 

The presence of a constant term in Table 5 is due to the invariance of the angular 
momentum. It would not be present in the reference frame of the invariant plane, 
orthogonal to the angular momentum of the system (see Malige et al, 2002). 
The existence of short periods is thus not an obstacle, and the aliasing prob- 
lem can be overcome, as in practice, using two output time steps is not a big 
constraint. It should not be costly in term of CPU time, as the two output are 
generated during the same integration. The large advantage of this procedure 
versus the numerical averaging that was usually performed online (as in No- 
bili et al, 1989), is that no information is lost, and the whole solution can be 
recovered. 
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